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ABSTRACT 

Detailed modeling of the millisecond brightness oscillations during thermonu- 
clear bursts from low mass X-ray binaries can provide important information 
about neutron star structure. Until now the implementation of this idea has 
not been entirely successful, largely because of the negligible harmonic content 
in burst oscillation lightcurves. However, the recent discovery of non-sinusoidal 
burst oscillation lightcurves from the accreting millisecond pulsar XTE J1814- 
338 has changed this situation. We, therefore, for the first time, make use of this 
opportunity to constrain neutron star parameters. In our detailed study of the 
lightcurves of 22 bursts, we fit the burst oscillation lightcurves with fully general 
relativistic models that include light-bending and frame- dragging for lightcurve 
calculation, and compute numerically the structure of neutron stars using realis- 
tic equations of state. We find that for our model and parameter grid values, at 
the 90% confidence level, Rc 2 /GM > 4.2 for the neutron star in XTE J1814-338. 
We also find that the photons from the thermonuclear flash come out through the 
layers of accreted matter under conditions consistent with Thomson scattering, 
and show that the secondary companion is a hydrogen burning main sequence 
star, with possible bloating (probably due to X-ray heating). 

Subject headings: equation of state — radiative transfer — relativity — stars: 
neutron — X-rays: binaries — X-rays: bursts 
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1. Introduction 

Nearly coherent brightness oscillations have been discovered with the Rossi X-ray Tim- 
ing Explorer (RXTE) during thermonuclear X-ray bursts from more than a dozen neutron 
stars in low-mass X-ray binaries (Strohmayer & Bildsten 2003). The large modulation am- 
plitudes at the onset of bursts, the time evolution of the pulsed amplitude during the rise of 
bursts (Strohmayer, Zhang, & Swank 1997), the coherence of the oscillations (Smith, Mor- 
gan, & Bradt 1997; Strohmayer & Markwardt 1999; Muno et al. 2000), and the long term 
stability of the oscillation frequencies (Strohmayer et al. 1998) led quickly to an interpreta- 
tion in terms of a "hot spot" on the surface of the star, so that the observed flux is modulated 
by rotation (Strohmayer et al. 1996). This interpretation has been strongly supported by 
observations of the transient accretion-powered pulsars SAX J1808-3658 (Chakrabarty et 
al. 2003) and XTE J1814-338 (Markwardt, Strohmayer, & Swank 2003), in which the fre- 
quency of burst oscillations is indeed extremely close to the spin frequency as inferred from 
persistent oscillations. 

Early on it was suggested that the intensity profiles of the pulse lightcurves contain 
valuable information about the stellar mass and radius (and hence about the equation of 
state of the dense matter in the stellar core), as well as the surface emission properties and 
system geometry (Strohmayer, Zhang, & Swank 1997; Miller & Lamb 1998; Weinberg, Miller, 
& Lamb 2001). The reason is that the lightcurves are affected by general relativistic light 
deflection as well as special relativistic beaming and aberration, which can be significant at 
the ~ 0.1 — 0.2c linear rotational speeds implied for spin frequencies of a few hundred Hertz. 
However, much of the required information is encoded in the ratios of pulse amplitudes at 
different harmonics of the spin frequency (although some constraints can be obtained if the 
amplitude at the fundamental is particularly high; see Nath, Strohmayer, & Swank 2002). 
Until recently, no source had harmonics detected definitively (for upper limits see Muno, 
Ozel, & Chakrabarty 2002). 

This situation has now changed. Analysis of XTE J1814-338, the fifth known accreting 
millisecond pulsar, shows a clear overtone in the pulse profiles of many individual bursts, 
with an amplitude that can be more than 0.25 times that of the fundamental (Strohmayer 
et al. 2003). In addition, the 314 Hz frequency of the fundamental oscillation is consistent 
with the spin frequency seen from the persistent pulsations. The existence of harmonics 
and the confirmation of the basic hot spot picture make this source very promising for de- 
tailed analysis and constraints on stellar, geometrical, and emission properties. We note that 
even more precise information is in principle available from the extremely accurately char- 
acterized lightcurves of the persistent pulsations of millisecond pulsars, which have several 
harmonics (e.g., see Cui, Morgan, & Titarchuk 1998; Poutanen & Gierlinski 2003). However, 
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the spectrum and profiles of pulses during accretion are affected by complicated details of 
shock physics and Comptonization. In contrast, thermonuclear burst spectra are well fit by 
blackbodies, without any trace of a high-energy tail related to Comptonization. They are 
also consistent with coming from very near the stellar surface. For these reasons, the pulses 
of XTE J1814-338 are most promising for an initial investigation. Combined with analy- 
ses of surface spectral line profiles (e.g., from EXO 0748-676; Cottam, Paerels, & Mendez 
2002), thermonuclear X-ray bursts are key sources for constraining the state of matter inside 
neutron stars. 

Here we report our analysis of the RXTE PCA data on 22 thermonuclear bursts from 
XTE J1814-338. In § 2 we briefly describe our data extraction techniques. In § 3, we discuss 
our theoretical models, the fitting procedure and the method of calculation of confidence 
intervals. We discuss our results in § 4. In particular, we show that for the two high- 
density equations of state (EOS) that we analyze, there is a preference for comparatively 
less compact neutron stars. We also show that the emission pattern from points on the surface 
is constrained moderately strongly and is close to that expected from the law of darkening 
for passive radiative transport in a Thomson scattering atmosphere (Chandrasekhar 1960). 
In this paper we use geometric units, i.e., G — c = 1. 



2. Data Analysis 

XTE J1814-338 was discovered in the Galactic bulge monitoring campaign carried out 
with the RXTE PCA (Markwardt & Swank 2003). The system has a binary orbital period 
of 4.28 hours and the neutron star has a spin frequency of 314.36 Hz (Markwardt et al. 2004, 
in preparation). This is the widest binary among the five known accreting millisecond pulsar 
systems. A total of 27 X-ray bursts were observed during the extensive RXTE follow-up 
of the outburst. A summary of basic burst properties can be found in Strohmayer et al. 
(2003). For the purposes of our modeling here we wished to produce lightcurves with as high 
a statistical precision as possible, while still maintaining energy resolved profiles. To do this 
we elected to phase-align different bursts (while preserving any misalignments among the 
energy bands). The energy channel boundaries were chosen so that the lightcurves in each 
band would have roughly similar statistical qualities (ie. similar numbers of total counts). 
We used the channel boundaries; 0-3, 4-6, 7-10, 11-13 and 14-28 of the "E_125us_64M_0_ls" 
PCA event mode. This corresponds to energy boundaries of ~ 2 - 3.7 keV, 3.8 - 5.0 keV, 5.2 
- 6.6 keV, 6.8 - 9.2 keV, and 9.4 - 23 keV. To co-add the bursts we found an offset phase for 
each burst which maximized the Z\ signal when adding the burst in question to the initial, 
reference burst. This is basically the same procedure described by Strohmayer & Markwardt 
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(1999). 

Although the different burst profiles are qualitatively similar, and hence adding them 
together is reasonable (see Figure 1 for added burst profile), we found that there is a general 
trend for the harmonic strength to decrease with time. That is, bursts occurring later in 
the outburst seem to have somewhat smaller ratios of the amplitude of the first overtone 
compared to that of the fundamental. This suggests that some secular change associated 
with the accretion rate might be influencing the details of the pulse profiles. Therefore, to 
allow for such changes while still obtaining good statistics, we added together the bursts 
in three groups based on their times of occurrence (see Figures 2, 3 and 4 for added burst 
profiles for three groups). Since the overall structure (ie. the mass and radius) of the neutron 
star must remain fixed through the outburst, the secular changes in the pulse profile may 
reflect small changes in the spot size, location or perhaps beaming function, with time. We 
will say more about this shortly. 



3. Physical Effects and Calculations 

The basic picture we employ is one in which after leaving the photosphere, photons 
propagate freely in vacuum to the observer. This approach ignores any scattering during the 
propagation. If there were scattering by a hot plasma near the star, we would expect this to 
leave its Comptonization imprint on the spectrum as an extended power-law tail. The lack 
of such a signature (Strohmayer et al. 2003) supports the assumption of free propagation. 

Even with this simplification, if the emission pattern on the surface and the angle- 
dependence of specific intensity from a given point on the surface are arbitrary, then there 
are too many free parameters and no meaningful constraints are possible. Therefore, to make 
progress in modeling, we adopt the following assumptions: 

• The surface consists of a background of uniform intensity and spectrum, plus a single 
hot spot that is a filled circle that emits uniformly. This is the most popular simple 
emitting system for burst oscillations. In a future paper we will consider models with 
two hot spots. 

• For R/M < 3.52 and for Schwarzschild spacetime, an emitted photon is deflected by 
more than 180° (Pechenick, Ftaclas, & Cohen 1983). Therefore, to simplify numerical 
techniques by ensuring that no emitted photon is deflected by more than 180°, the 
minimum value of R/M we consider is 3.6. 

With these assumptions, we calculate lightcurves by tracing rays backwards, i.e., from 
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the observer to the stellar surface. The paths of these rays, and their specific intensities, carry 
the imprint of several physical effects. Using an approach similar to that of Ozel & Psaltis 
(2003), we consider: (1) Doppler boosts, (2) special relativistic beaming, (3) gravitational 
redshift, (4) light-bending, and (5) frame-dragging. We do not include the effects of spin- 
induced stellar oblateness. These effects are only second order in the stellar rotation and 
are thus small compared to other uncertainties. For example, for all our EOS models, the 
polar to equatorial radii ratio of the star is greater than 0.98 (for stellar mass M > 1.4M 
and observed stellar spin frequency ~ 314 Hz). We describe tests of our raytracing and 
lightcurve codes in the Appendix. 

To make our model lightcurves realistic, we calculate the relation between M and the 
stellar radius R, as well as a/M (where, a is the stellar angular momentum per unit mass) 
for a given M and (~ 314 Hz for XTE J1814-338), with a specified EOS model. We 
do this by computing numerically the stable structure of a rapidly spinning neutron star, 
using the formalism given by Cook, Shapiro, & Teukolsky (1994) (for a detailed description, 
see also Bhattacharyya et al. 2000; Bhattacharyya 2002). For a particular neutron star 
EOS model and assumed values for the stellar mass and spin rate, we solve Einstein's field 
equations and the equation of hydrostatic equilibrium self-consistently to obtain other stellar 
structure parameters (radius, angular momentum, etc.). We use these output parameters in 
our timing calculations. 

In our detailed lightcurve calculation, we have five input parameters, for a given EOS 
model and the known value of v*. In our study, the structure of the star is fixed by one 
unknown parameter, which we choose to be (1) the dimensionless stellar radius to mass ratio 
(R/M). The hot spot is specified by two parameters: (2) the ^-position 9 C of the center of 
the spot, and (3) the angular radius A6 of the spot. The emission from a single point on the 
spot, as measured in the corotating frame, is characterized by (4) a parameter n that gives a 
measure of the beaming in the emitter's frame (corotating with the star), where the specific 
intensity as a function of the angle ip (in the emitter's frame) from the surface normal is 

oc cos n (ip). Finally, we have (5) the inclination angle i of the observer as measured from 
the upper rotational pole. It is to be noted, that at the initial phase of our model calculation, 
we consider that the hot spot emits as a blackbody (as burst spectra can be fitted well by 
blackbodies) with a temperature kT = 2 keV. However, later we allow the source spectrum 
to deviate from a blackbody (for discussion about the spectrum and the justification of the 
assumed temperature, see later in this section). 

In our work, we consider two representative EOS models: A18 and Al8+<5f+UrX (Ak- 
mal, Pandharipande, & Ravenhall 1998). The first one is the Argonne t>i 8 model (Wiringa, 
Stoks, & Schiavilla 1995) of two-nucleon-interaction. This model fits the Nijmegen database 
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(Stoks et al. 1993) with xV^data ~ 1, and hence is called "modern" . The model Al8+fo;+UIX 
considers two additional physical effects: the three-nucleon-interaction (Urbana IX (UIX) 
model; Pudliner et al. 1995) and the effect of relativistic boost corrections to the A18 in- 
teraction. The EOS model A18 is soft, i.e., the maximum non-rotating mass it supports is 
small (~ 1.67M ). On the other hand, the EOS model A18+5f+UIX is hard (maximum 
supported non-rotating mass is ~ 2.2M Q ). Figure 2 shows the stable stellar mass vs. radius 
curves for these EOS models (for stellar spin frequency ~ 314 Hz). Although, there are 
many other EOS models available in the literature, our chosen models span a representative 
range. 

In order to compute the energy dependent flux, we trace back the paths of the photons 
from the observer to the hot spot, using the Kerr spacetime. We solve the geodesic equations 
numerically using a fifth order Runge-Kutta method with adaptive stepsize control. The 
accuracy of the results is tested and described in the appendix. Our ray-tracing method is 
similar to that of Bhattacharyya, Bhattacharya, & Thampan (2001), except that here we use 
the Kerr spacetime, instead of the correct spacetime around a rapidly rotating neutron star, 
for the sake of simplicity. It is unlikely that the errors introduced by this approximation 
will have any detectable effect on the lightcurves for the stellar spin frequency we consider 
(i.e., ~ 314 Hz). Our method has two major differences from the approach of, e.g., Braje, 
Romani, & Rauch (2000). (1) We solve the geodesic equations of motion of a photon from 
the observer to the surface of the neutron star, instead of from the surface to the observer. 
(2) Instead of doing elliptic integral reduction of the geodesies (as was done by Braje et al. 
2000), we solve the equations of motion numerically, for easier generalization to numerical 
spacetimes. 

The ray-tracing enables us to get the boundary of the image of the source at the ob- 
server's sky. We then calculate the observed flux by integrating the observed energy de- 
pendent specific intensity inside the boundary of this image (for a detailed description, see 
Bhattacharyya et al. 2001). The model lightcurve is calculated by repeating the same pro- 
cedure for many spots at different 0-positions (but the same ^-position) on the surface of 
the star. The actual phase points of the lightcurve are calculated from these 0-positions, the 
stellar spin frequency, and the time delay consideration. The time delay is because photons 
emitted at different points on the stellar surface take different times to reach the observer. 

In this paper, we compare our models with the data of 22 bursts. Based on burst 
occurrence time and the harmonic strength of the lightcurves, as described above, we added 
all the bursts in three groups; bursts 1-8; 9-16; and 17-22. We phase-align and stack the 
bursts within each group, and therefore effectively analyze three characteristic sets of pulse 
profiles. We used five energy bands for each group. The energy ranges are given in §2 above. 
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We compare our models directly with the data. To do this, we compute a model 
lightcurve of intensity and spectrum versus pulse phase. In this initial step we consider 
only the lightcurve of the hot spot, not any background emission. We fold the model spec- 
trum through the RXTE response matrix to get counts per channel as a function of rotational 
phase. For each channel range, we then add a phase-independent background (which is a 
free parameter) and normalize it so that the total number of counts in that channel range, 
summed over all phases, is the same as the number of counts in that channel range in the 
observed lightcurve. Our final step is to shift the entire lightcurve in phase (by same amount 
for all the channel bands) to maximize the match with the observed lightcurve, as determined 
by a x 2 statistic. 

The procedure of adding backgrounds and normalizing counts independently in each of 
the channel ranges has two implications. First, the independent backgrounds mean that we 
leave the unpulsed spectrum unconstrained. Second, the independent normalizations mean 
that although we calculate the initial spectrum and lightcurve of the hot spot assuming 
a kT = 2 keV blackbody, renormalization allows the effective spectrum to deviate from 
blackbody. Therefore, we preserve the rms variability of a blackbody with, e.g., Doppler 
shifts, but not the whole spectrum. Indeed, in reality we do not expect the temperature to 
be constant throughout a hot spot; for example, if the hot spot is related to a magnetic pole, 
it might be hotter in the center than at the edges. In addition, the pulse profiles extracted 
from the observations incorporate data over a period of several tens of seconds, during which 
the hot spot is likely to cool significantly. Thus, the time-averaged spectrum will not truly be 
a blackbody, hence we let the spectrum be adjusted by allowing independent normalizations. 

The initially chosen hot spot temperature of 2 keV is justified for the following reasons: 
(1) It is expected to be close to the average blackbody temperature of the hotspot for most 
of the bursts (see Kuulkers et al. 2003 for blackbody temperature variation during burst 
evolution for several sources), and (2) The likelihood distributions (see next two paragraphs 
for description) for the parameters are similar, even for widely different initially chosen hot 
spot temperatures. 

Using the above procedure of adding backgrounds, normalizing, and shifting phases, we 
find the best \ 2 value for many different combinations of parameters. To be systematic, 
we consider combinations of parameters on a grid (for the values of each parameter, see 
Table 1). We also use the Markov Chain Monte Carlo (MCMC) method to ensure that 
the grid method does not miss any significantly low x 2 values. The number of counts per 
phase bin in each channel range is large enough that Gaussian statistics apply, hence the 
likelihood is proportional to exp[— x 2 /2]. We apply the physical constraint that the value of 
each parameter is same for all the channel ranges for a given burst group (some parameters, 
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for example, the size and the location of the hot spot may change from one burst group 
to another). Therefore, for each parameter combination, the likelihood is proportional to 
exp[— S% 2 /2], where the sum is over all the channel ranges. 

Once the likelihoods are computed for each combination of parameters, we produce 
likelihoods for each parameter individually via the process of marginalization. In this pro- 
cess, let the posterior probability density over the full set of model parameters 9± . . .0 n be 
p(6i . . . 6 n ). If we are interested in the credible region for a single parameter 6^, then we inte- 
grate, or marginalize, this probability distribution over the "nuisance" parameters 9i . . . Qk-\ 
and 9 k+1 . . . 6 n : 



For our purposes we assume that each combination of parameters has the same a priori 
probability, hence the posterior probability is simply proportional to the likelihood. For 
Gaussian distribution of likelihoods of a parameter, this method is equivalent to the standard 
frequentist procedure. However, for non-Gaussian distribution, this Bayesian method is more 
general. 

We find that the posterior probability distributions are not Gaussian for most of our 
parameters. Therefore, although there is a unique maximum likelihood value for any param- 
eter, credible regions cannot be quoted uniquely. For example, one gets a different value if 
one assumes a symmetric distribution than if one looks for the smallest region containing a 
certain total probability. We have adopted a variant of the latter procedure. Using linear 
interpolations of the probability between grid points for a given parameter, we compute the 
smallest region that includes the maximum likelihood point and contains 90% of the to- 
tal probability. For parameters expected to remain unchanged between bursts (specifically, 
R/M and the observer inclination i, and also the beaming parameter n) we combine data 
from all bursts. For the other parameters, we calculate likelihoods separately for each of the 
three burst groups. 



As discussed in § 3, we perform comparisons using two EOS models. Figure 6 shows 
the goodness of fit for a certain burst group and channel range, and for a representative 
combination of parameter values. The x 2 value for this example comes out to be ~ 25.5, for 
the number of degrees of freedom equal to 24. This goodness of fit suggests that the overall 
model framework is a reasonable representation of the data. 

In Figure 7, we plot the likelihoods with the parameter R/M, for two EOS models, 




(1) 



4. 



Results and Conclusions 



- 9- 



using the data from all 22 bursts. While we can not determine an upper limit of R/M, lower 
limits (4.7 for A18+Sv +UIX and 4.2 for A18) of 90% confidence regions (given in Table 2) 
can be computed. We do not consider R/M values greater than 6.0, as for these values, the 
stellar mass is too small (< l.O7M for A18 and < 1.29M for Al8+5w+UIX) to be realistic 
for an accreting neutron star. A simple extrapolation of the likelihood curves in Figure 7 
shows that the probability of R/M < 3.6 is very small. However, for R/M < 3.52 and for 
Schwarzschild spacetime, photons from a single point on the stellar surface may reach the 
observer by multiple paths. As a result, the corresponding lightcurves may be qualitatively 
different from the ones we calculate. Therefore, we can say that one interpretation of the 
data is that R/M is greater than 4.2 with 90% confidence, but we can not completely exclude 
the possibility of R/M < 3.52. 

Figure 8 displays the likelihood distribution with the observer's inclination angle i. 
This figure shows that i > 22° with very high probability. Here, as well as for the likelihood 
calculation of other parameters, we focus our fitting procedure on the range 20° — 50° for i. 
This is because, with a smaller number of grid points for other parameters, our fitting for 
inclination angle values in the range 5° — 90° shows that the likelihood values for % < 20° 
are very small, and the inclusion of inclination angles greater than 50° produces only minor 
changes to our results. However, constraints on inclination angle by observations in other 
wavelengths would help us constrain other parameters more effectively. Presently, M. Krauss 
et al. (2004, in preparation) are working on this using optical data. 

We keep the beaming parameter n (defined in § 3) fixed in all three burst groups in 
our statistical procedure. This is because the thermonuclear flash is expected to occur at 
the bottom of a pile of accreted matter. Therefore, although the original emission may be 
isotropic, what comes out through the layers of accreted matter should be beamed. As the 
optical depth of these layers is expected to be very high, we expect the law of darkening 
(and hence the value of n) to be the same for all the bursts. 

The likelihood distributions for n for two EOS models are shown in Figure 9. For both 
EOS models, the beaming parameter n peaks at n — 0.75. This figure shows that n is well 
constrained on lower side, and the lower limit of the 90% confidence regions (for both the 
EOS model) is 0.55 (given in Table 2). The law of darkening for semi-infinite plane-parallel 
layers with a constant net flux and Thomson scattering is given in Chandrasekhar (1960). 
In Figure 10, we compare this law of darkening with our emission function cos n (ip). 
We find that, for n = 0.55 (where likelihood value is ~ 15 — 20% of the peak likelihood 
value), our emission function matches with I{jp) quite well, except for nearly tangential 
emission (for which the emitted flux is small). The matching is also reasonably good for the 
most probable value n = 0.75. Therefore, the data are consistent with Thomson scattering 



-10- 



through an optically thick layer (probably in the accumulated accreted matter) of thermal 
electrons. 

In Table 3, we show the 90% confidence intervals for the polar angle of the center of the 
hot spot (9 C ) for both the EOS models. Here we give the results for the three burst groups 
separately, as the values of this parameter may change from one group to another. The union 
of 90% confidence intervals for all burst groups and EOS models is 60° — 139°. Assuming 
that the hot spot is at the magnetic pole, this implies an angle of 40° — 90° between the 
spin axis and the magnetic axis. The peak of the likelihood distribution for 6 C always occurs 
around the interval 90° — 110° (see Figures 11, 12 & 13), hence the magnetic pole may be 
close to the rotational equator. 

Table 4 displays the likelihood values for the angular radius (A9) of the hot spot for 
different burst groups and EOS models. For the burst group 1 — 8, A9 is comparatively 
better constrained, and a smaller spot (A6 ~ 5° — 25°) is likely. For the other two burst 
groups, comparatively larger spots are probable. 

A direct comparison of the likelihoods of two EOS models shows that the difference is 
not enough to constrain EOS models by this direct method. 

The pulsar mass function and orbital period for this source are 0.002016 M Q and 4.27 
hours respectively (Markwardt, Strohmayer, & Swank 2003). Hence, a reasonable range of 
neutron star mass (1.4 M — 2.0 M ) and our % > 20° inclination angle constraint imply 
a companion mass of 0.17 M — 0.72 M & . For this range of masses and for the observed 
orbital period, the companion can be neither a helium main sequence star nor a degenerate 
star (Bhattacharya & van den Heuvel 1991). It is also too large to be a brown dwarf (as 
was predicted for another accreting millisecond pulsar system SAX J1808. 4-3658; Bildsten & 
Chakrabarty 2001). The most probable option is that it is either a hydrogen burning main 
sequence star, or an evolved (sub-)giant star. However, the companion mass is too small for 
it to have evolved off the main sequence (Bohm-Vitense 1992). Therefore, the most likely 
option is a hydrogen main sequence star. This is also seen in Figure 14. If the companion 
is a normal hydrogen burning main sequence star (as shown in Figure 14), its maximum 
possible mass is ~ O.48M . This is because, the radius of the companion can not be smaller 
than that of a hydrogen main sequence star (for a given mass), but it can be bigger, because 
it can be bloated due to X-ray heating by the neutron star and/or the accretion disk. This 
upper limit of companion mass corresponds to a lower limit of observer's inclination angle 
24° (using the pulsar mass function value and a lower limit of neutron star mass 1.4M ). 
Interestingly, this lower limit of inclination angle is close to that found by our lightcurve 
fitting. However, if the mass fraction of hydrogen in the hydrogen burning main sequence 
companion star is very low, then it may have a higher mass than O.48M (see Rappaport & 
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Joss 1984). If we consider the lower limit (0.17 M ) of the companion mass, then its radius 
is about 98% bigger than the normal radius (Bhattacharya & van den Heuvel 1991). 

In conclusion, the significant overtones observed in the pulse profiles of burst brightness 
oscillations from XTE J1814-338 open up new possibilities for constraints on neutron star 
structure as well as on source geometry and emission properties. In a subsequent paper, we 
will study the analysis possible with future large-area detectors, in particular whether EOS 
models may be constrained strongly based on their likelihoods. 
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We use three codes in series to calculate the \ 2 values. These codes have been checked in 
different ways. Here we mention a few of them. The ray-tracing code traces back the paths 
of photons from the observer to the surface of the star, and calculates different parameter 
values at the footprints of the photons on the surface. For a photon that is emitted from the 
surface tangentially (i.e., the one with maximum angular momentum), one can calculate the 
amount of angular deflection at infinity (for the Schwarzschild spacetime). The comparison 
between these values with the outputs from our code give satisfactory results. For example, 
when R/M = 4.2, the total angle traveled by a photon emitted tangentially at the surface is 
145°. 8, and our code yields 145°. 1. The code also calculates the total redshift with typically 
less than 0.01% error. However, a photon travelling directly over a pole of the star gives a 
few degrees of error in the ^-position of its emission point. But the error in this single ray 
does not introduce significant error in our final results. 

For both the Schwarzschild and Kerr metrics, the value of L = —p^/pt can be derived 
analytically at the observer's position and at the surface of the star. Here p<f> and — p t are 
photon's 0-angular momentum and energy respectively. At the observer's position 



where, b and a are plane polar coordinates in observer's sky, and i is the observer's inclination 
angle. At the surface of the star, 



Appendix 



L = — b sin i sin a , 



(2) 



L = 



9<w cos£ 



(3) 




This preprint was prepared with the AAS IATgX macros v5.2. 
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where, £ is the emission angle with ^-direction, and g^'s are the metric coefficients of a 
metric 

dS 2 = g tt dt 2 + g„dr 2 + g ee dd 2 + g H d<\> 2 + g^d^dt . (4) 

As L remains constant along the photon's path, the values of L from equations (2) and (3) 
should be the same. We find that with our ray-tracing code the maximum difference in these 
two values is ~ 1%. 

The lightcurve calculation code computes the energy dependent observed flux from 
different hotspots on the stellar surface. This is done by the integration of observed specific 
intensity over the image of the hotspot at the observer's sky. As the image of the hotspot 
is of irregular shape, this integration is done using a Monte Carlo method. We ensure 
that different random seed values and different number of points (used for Monte Carlo 
integration) give a stable result. For example, for parameter combinations giving x 2 /dof ~ 1, 
differences in x 2 s are less than 0.5%. For larger x 2 values, these differences are more, but 
still very low, and these have little effect because it is the lower x 2 values that dominate the 
likelihoods. These differences in x 2 s appear to be more or less random, hence the differences 
in total x 2 s (addition for different channel ranges and burst groups) are also small. The 
integrated area is also verified to be correct (with typically less than 0.1% error) even for a 
small spot (for example, with 5° angular radius) with its center at the same and 6 values 
as of the observer. 

Finally, the comparison code compares the theoretical lightcurves with the observed 
data, and computes the x 2 values. We check this code in different ways. For example, we 
consider different maximum background values and different steps in background values to 
ensure a stable result (maximum ~ 0.03% difference in x 2 values). We also use synthetic 
data to confirm that the known parameters are obtained to accuracies compatible with our 
uncertainty region for the real data. 
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Fig. 1. — The fully added observational lightcurve (22 bursts and — 28 channels). 
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Fig. 2. — The observational lightcurve (bursts: 1 — 8 and — 28 channels). 



-17- 



73 

o 



2.20xl0 4 
2.10xl0 4 

2.00xl0 4 

1.90xl0 4 



u l.BOxlO 4 



1.70xl0 4 



1.60xl0 4 



l.SOxlO 4 




0.0 



0.2 0.4 0.6 0.8 

Phase (arbitrary) 



1.0 



Fig. 3. — The observational lightcurve (bursts: 9 — 16 and — 28 channels). 
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Fig. 4. — The observational lightcurve (bursts: 17 — 22 and — 28 channels). 
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Fig. 5. — Mass vs. radius diagrams: Curve 1 is for the EOS model Al8+5v+UIX and curve 
2 is for the EOS model A18. These diagrams are for stable stellar configurations, spinning 
with a frequency ~ 314 Hz. The solid parts of the curves (correspond to higher R/M sides of 
the vertical lines in Figure 7) are the allowed regions with 90% confidence interval of R/M. 
These regions constrain the mass vs. radius plane for the neutron star in XTE J1814-338. 
The solid oblique line corresponds to R/M = 3.6. In this work, our results are for the M — R 
space below this line. 
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Fig. 6. — An example fit to the data. Stacked data (solid curve) for bursts 9 — 16 (channel 
range: 7 — 10) from XTE J1814-338, versus a model (dotted curve) in which R/M = 4.9, 
i = 36°, 9 C = 110°, the hot spot has an angular radius of 45°, and n = 0.8 (see text for 
description of parameters). Here we have used the EOS model Al8+5f+UIX. In this fitting 
we get a x 2 value ~ 25.5, for the number of degrees of freedom equal to 24. 



- 21 - 




Fig. 7. — Likelihood distribution of stellar radius to mass ratio R/M, using data from all 22 
bursts. The solid curve is for the EOS model Al8+5v +UIX, and the dashed curve is for the 
EOS model A18. The solid vertical line gives the lower limit of the 90% confidence region 
for the EOS model Al8+Sv +UIX, and the dashed vertical line gives that for A18. 
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Fig. 8. — Likelihood distribution of observer's inclination angle i, using data from all 22 
bursts. Solid curve is for the EOS model Al8+5v +UIX, and the dashed curve is for the 
EOS model A18. The solid vertical lines give 90% confidence interval for the EOS model 
Al8+<k> +UIX, while the dashed vertical line gives the lower limit of 90% confidence region 
for A18. This figure demonstrates that a value of % less than 22° is highly improbable. 
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Fig. 9. — Likelihood distribution of the beaming parameter n (defined in § 3), using data 
from all 22 bursts. The solid curve is for the EOS model Al8+5t>+UiX, and the dashed 
curve is for the EOS model A18. The solid vertical lines give 90% confidence interval for 
the EOS model Al8+<k> +UIX, while the dashed vertical lines give that for A18. This figure 
demonstrates that a value of n less than 0.2 is highly improbable. 
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Fig. 10. — Comparison between (Law of darkening; Chandrasekhar 1960) and our emis- 
sion function (cos™ ip) : Both the functions give the angular distribution of specific intensity 
in the emitter's frame (corotating with the star). Here the independent parameter is the 
emission angle (tp) (measured from surface normal) in the emitter's frame, and the depen- 
dent parameter is the ratio of Acos n (-?/>) to I(ip) (where, A is an arbitrary constant, and I(ip) 
is the law of darkening for the radiative transfer in a semi-infinite plane-parallel atmosphere 
with a constant net flux and Thomson scattering). The solid curve is for n = 0.55, and the 
dashed curve is for n = 0.75. For the former one, is well fitted by our emission function, 
except for very high emission angle, for which the emitted flux is small. 
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Fig. 11. — Likelihood distribution of the polar angle of the center of the spot (9 C ) for bursts 
1—8. Solid curve is for the EOS model Al8+5v +UIX, and the dashed curve is for the 
EOS model A18. The solid vertical lines give 90% confidence interval for the EOS model 
AlS+Sv +UIX, while the dashed vertical lines give that for A18. This figure demonstrates 
that a value of 9 C less than 50° or greater than 150° is highly improbable for bursts 1 — 8. 
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Fig. 12. — Likelihood distribution of the polar angle of the center of the spot (9 C ) for bursts 
9 — 16. Meanings of all the curves and lines are same as in Figure 11. This figure demonstrates 
that a value of 9 C less than 50° is highly improbable for bursts 9 — 16. 
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Fig. 13. — Likelihood distribution of the polar angle of the center of the spot (6 C ) for 
bursts 17 — 22. Meanings of all the curves and lines are same as in Figure 11. This figure 
demonstrates that a value of 6 C less than 50° is highly improbable for bursts 17 — 22. 
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Fig. 14. — Companion mass vs. radius plane (similar to the Figure 4 of Markwardt et al. 
2002): The lines under the phrase 'XTE J1814-338' are for the companion star in the source 
XTE J1814-338 (for an observer inclination angle > 20°); the dashed line is for a neutron 
star mass 1.4M and the dotted line is for the neutron star mass 2.OM . For the purpose of 
comparison, we plot the mass-radius curves for hydrogen burning main sequence star (solid 
line; using R/R Q = M/M Q ), helium main sequence star (dash-dot line; Verbunt 1990), 
brown dwarfs (dashed lines), and cold helium dwarf (dotted line). Brown dwarf models are 
for ages of 0.1, 0.5, 1, 5, and 10 billion years (from top to bottom). This figure demonstrates 
that the companion star in the source XTE J1814-338 is likely to be a hydrogen burning 
main sequence star. 
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Table 1. Grid values considered for the parameters for likelihood calculations. 



Parameter 


Grid values 




3.60, 3.70, 3.80, 3.90, 4.00, 4.10, 




4.15, 4.20, 4.25, 4.30, 4.35, 4.40, 


R/M 


4.45, 4.50, 4.55, 4.60, 4.65, 4.70, 




4.75, 4.80, 4.85, 4.90, 5.00, 5.10, 




5.20, 5.30, 5.40, 5.50, 5.60, 5.70, 




5.80, 5.90, 6.00 




20°, 22°, 24°, 26°, 28°, 30°, 32°, 


i 


34°, 36°, 38°, 40°, 42°, 44°, 46°, 




48°, 50° 




10°, 30°, 50°, 70°, 90°, 




110°, 130°, 150°, 170° 


A9 


5°, 25°, 45°, 60° 




0.00, 0.20, 0.40, 0.45, 0.50, 


n 


0.55, 0.60, 0.65, 0.70, 0.75, 




0.80, 0.90, 1.00, 1.10, 1.50 
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Table 2. 90% confidence intervals for R/M, i and n. 



Parameter 


EOS 


Lower limit 


Upper limit 


R/M 


A18+5v+XJIX 


4.7 


6.0 




A18 


4.2 


6.0 


i 


A18+5v+VlX 


26° 


48° 




A18 


36° 


50° 


n 


A18+5v+VlX 


0.55 


1.30 




A18 


0.55 


1.17 
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Table 3. 90% confidence intervals for the polar angle (9 C ) of the center of the spot. 



Parameter EOS Burst group Lower limit a Upper limit a 







Bursts: 


1 - 


8 


69° 


131° 




A18+<to+UIX 


Bursts: 


9- 


16 


82° 


138° 






Bursts: 


17- 


- 22 


82° 


139° 






Bursts: 


1 - 


8 


60° 


118° 




A18 


Bursts: 


9 - 


16 


74° 


129° 






Bursts: 


17- 


- 22 


69° 


129° 



a We calculate the 90% confidence interval for each burst group separately, 
because this parameter could in principle vary between bursts. 
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Table 4. Likelihood distribution of the angular radius of the spot AO. 



Burst group 3. 


Parameter 


EOS: A18+5v+VlX 


EOS: A18 




AO 


Likelihood 


Likelihood 




5° 


1.000 


1.000 


Bursts: 1 — 8 


25° 


0.566 


0.426 




45° 


0.503 


0.082 




60° 


0.047 


0.002 




5° 


0.718 


0.826 


Bursts: 9-16 


25° 


0.458 


0.712 




45° 


1.000 


1.000 




60° 


0.379 


0.153 




5° 


0.367 


0.569 


Bursts: 17-22 


25° 


0.439 


0.649 




45° 


1.000 


1.000 




60° 


0.354 


0.110 



a We calculate likelihood distribution for each burst group sep- 
arately, because this parameter could in principle vary between 
bursts. 



